Computer Algebra Math0070

Daniel Richardson, Department of Mathematical Sciences, University of Bath email:dsr@maths.bath.ac.uk

May 15, 2001

1 Introduction: Computer Algebra

To see the latest version of these notes, from BUCS, in unix command:

xdvi ~masdr/public_html/m70

Computer algebra is the use of computers to represent mathematical objects and to compute with them. See Computer Algebra, J Davenport for discussion. The mathematical objects include numbers, but also polynomials, sets, differential equations , etc. Examples of computer algebra systems are maple, reduce, axiom. We will start with maple.

1.1 Maple

There are a number of books in the library which are intended to help with use of Maple. See (518.542MAP). From BUCS network in unix, command xmaple and then choose help and do the introductory tour. Try some commands.

Example 1 You might try the following. (Comments in parentheses.)

1/3+2/7; (calculator) a:= sqrt(2); ( a now stands for the real number sqrt(2), not an approximation) a; (what do you expect) evalf(a); (evalf gives floating point approx with current precision, held as Digits) Digits:=50; ( note that Digits and digits are not the same) evalf(a); (50 digits precision) evalf(Pi); (note that Pi and pi are not the same) p:= x+y; (p now stands for x+y) p; q:=p**5; (You can also use carat for exponentiation) q:= x*q + y*q; expand(q); (multiply and collect) q/p**10; simplify(%); (percent sign means pervious value, simplify tries to simplify) subs(x=4.5, q); ( substitution) sort(q, [x,y]); sort(q, [y,x]); with(linalg); (turns on the linear algebra package) b:=matrix([[x,y], [x*z, z]]); b; evalm(b); ( displays matrix b) c:= 1/b; (matrix inverse) d:= b&* c; ( &*, ampersand followed by *, is matrix multiplication) d; simplify(d); simplify(evalm(d)); (only at this stage do we get identity displayed) subs(x= 2.1, evalm(c)); ( matrix substitution, note evalm inside subs) l1:=[a,b,d]; ( list, refer to help) s1 := {a, b, c}; ( set, order independent, no duplicates, as mathematical set) sq:= a,b,c; ( sequence) plot(sin(x)*exp(-x/5), x= 0 .. 4*Pi); (try also 3d plots)

Exercise 1 This is an exercise in interactive use of maple. It is harder than it looks. Please persist. The problem is to tell maple what we want it to do. Note that the whole computation is somewhat ambiguous, and there is room for doubt about the correctness of the computation. What is the correct framework? What are x_{1}, x_{2}, x_{3}? Real? Which square root is meant? The only way to be sure that this result is correct is to make a proof with pencil and paper.

First turn on linear algebra package by with(linalg); Define x to be the 1 by 3 row matrix [x_{1}, x_{2}, x_{3}]. For j between 1 and 3, define s_{j} to be (∑_{j≤k≤3} x_{k}^{2})^{1/2}. Define H to be a 3 by 2 matrix, with H[i, j] = 0 if i < j H[i, i] = s_{i+1}/s_{i} for i = 1, 2 H[i, j] = −x_{i}x_{j}/(s_{j}s_{j+1} for 1 ≤ j < i ≤ 3. Verify:

1. transpose(H)H is 2 by 2 identity

2. xH is 0.

3. Let p = H ∗ transpose(H). Then transpose(p) = p, and p ∗ p = p

The matrices H and p are basic parts of the PSLQ algorithm, which will be discussed later.

1.2 Procedures in Maple

Advice: make a directory to hold all maple files and go into this directory before running maple. Create maple source with any editor, e.g. emacs. Save it, e.g. fb.mpl. Start maple. Command.

read “fb.mpl”; and then use procedure.

Example 2

fib := proc(x:: nonnegint)
# Fibonacci function
option remember;
if x = 0 then 0
elif x = 1 then 1
else fib(x-1) + fib(x-2) fi
end:

# some tests
fib(1);
fib(6);
fib(7);
fib(8);

These notes will gradually expand and will be handed out gradually. A solution to the first exercise is at

~masdr/public_html/exer1.mws

and you should be able to open it from Maple. For the problems class on Thursday, try the following:

Exercise 2 Consider the above Fibonacci procedure and the two following which were discussed in class.

cheb := proc(n:: integer)
# returns table of Chebyshev polynomials up to n
  local p,k;
  p[0]:=1;
  p[1]:=x;
  if n <= 1 then RETURN(eval(p))
  fi;
  for k from 2 to n do
    p[k] :=expand(2*x*p[k-1] - p[k-2])
  od;
  eval(p)
end:

lop := proc(f::algebraic,g::algebraic,x::name, a)
#  this procedure attempts to find the limit, as x-->a of
#  f(x)/g(x), using L'Hopital method.  It assumes f,g are expressions
#   for infinitely differentiable
#   real valued functions of a real variable called x
  if eval(subs(x=a,g)) <> 0 then eval(subs(x=a,f/g))
  elif eval(subs(x=a,f)) <> 0 then "no finite limit"
  else lop(diff(f,x), diff(g,x),x,a)
  fi
end:

1. Find limit_{x→0}(x^{3}sin(x)/(1 − cos(x))^{2}.

2. Modify the Fibonacci procedure to return a table.

3. Modify the Chebyshev procedure to return the nth polynomial

4. Modify the L’Hopital procedure so that it does not use recursion.

5. Does the L’Hopital procedure always terminate? Can you fix or alleviate this problem?

6. Consider the Euclid GCD algorithm given below. Write a procedure for the extended Euclidean algorithm, which returns a,b,g where g is the gcd(x,y) and g=ax+by.

2 Problems and Algorithms

Notation: N, Z, R, C are, respectively, the natural numbers including 0, the integers, the real numbers and the complex numbers. Z[x] is the set of polynomials in x with integral coefficients, i.e. expressions of the form a_{0}x^{d} + ... + a_{d}, where a_{0}, ...a_{d} are integers. Z[x_{1}, ..., x_{n}] is the set of polynomials in variables x_{1}, ..., x_{n} with integral coefficients.

Problem means a family of related problems. Examples:

1. GCD over Z. Given integral x, y, find their greatest common divisor.

2. Factorisation over Z. Given n in Z, express n as a product of prime numbers.

3. Factorisation over Z[x], down to integers. Given p(x) in Z[x], express p[x] as a product of integers and irreducible factors in Z[x].

4. (Hilbert’s tenth Problem). Given p(x_{1}, ..., x_{n}) in Z[x_{1}, ..., x_{n}], decide whether or not there exist integers z_{1}, ..., z_{n} so that p(z_{1}, ..., z_{n}) = 0.

5. (Polynomial System Solving) Given (p_{1}, ..., p_{n}) : C^{m} → C^{n}, and each p_{i} in Z[x_{1}, ..., x_{m}], decide whether or not there exists α in C^{m} so that (p_{1}, ..., p_{m})(α) = 0.

6. Identity problem for trigonometric and absolute value Functions. Given expression for function f(x) built up from x, π, 1, +, −, ∗, sin(x), cos(x), | x | and composition, decide whether or not f(x) is identically zero.

7. Integer Relation Finding. Given real numbers x_{1}, ..., x_{n} decide whether or not there exist integers m_{1}, ..., m_{n} not all zero and with absolute value bounded by M so that m_{1}x_{1} + ... + m_{n}x_{n} = 0.

An algorithm to solve such a problem is a method which given an instance of the problem will determine a solution after finitely many steps. Assume some measure of the size of the input. We measure the complexity of an algorithm as the number of steps it takes as a function of the size of the input. We will get different versions of complexity, depending on what we allow as a single step.

1. Classical bit complexity. Number of steps of a Turing machine.

2. Complexity over Z. We allow addition, subtraction and multiplications of arbitrary integers in one step.

3. Complexity over R. We allow addition, subtraction, multiplication, and testing of sign for arbitrary real numbers in one step. Some problems, e.g. Integer relation finding as above, only make sense over R.

We can classify problems by difficulty:

1. Recursively unsolvable. No algorithm is possible. Hilbert’s tenth problem and the Identity problem for trigonometric and absolute value functions are examples. There is no hope of a uniform solution method.

2. Recursively solvable, but intractable. In this case, an algorithm exists, but the complexity of the solution takes a number of steps which increases faster than any polynomial in the input size. Most people believe that factorisation over Z is in this category (with size of number measured by length of digital representation). Note that even for an intractable problem, a computer program may be useful.

3. Polynomial complexity solvable. These are problems for which an algorithm exists with the number of steps for a solution bounded by some polynomial in the size of the input. Integer relation finding is in this category, with the size of the bound M measured as the number of digits to express M, and with complexity over R. Factorisation over Z[x] down to integers is another problem in this category.

4. Feasible in practice. GCD over Z for example.

See www.computer.org/cise/cs2000/c1toc/htm

for a list of “the best algorithms of the century”. It is a good idea to read all these articles, but we are only going to study one of them.

If possible we would like to translate algorithms into usable computer programs. A computer algebra system such as Maple is a good framework for this. We will consider integer GCD, Integer relation finding, numerical solving of systems of equations by Newton’s method, and factorisation over Z[x].

3 Euclid’s GCD for Z

For integers x and y, x | y means that x is a factor of y. An integer d that divides both x and y is called a common divisor of x and y. GCD(x, y) is the greatest common divisor of x and y. To find GCD, we can form the remainder sequence:

n_{0} = x

n_{1} = y

n_{0} = q_{0}n_{1} + n_{2}, where 0 ≤ n_{2} <| n_{1} |

n_{1} = q_{1}n_{2} + n_{3}, where 0 ≤ n_{3} <| n_{2} |

...

n_{i} = q_{i}n_{i+1} + n_{i+2}, where 0 ≤ n_{i+2} <| n_{i+1} |

...

n_{k} = q_{k}n_{k+1} + 0

until n_{k+2} = 0. Then n_{k+1} is the GCD Of x and y. In each case, n_{i+2} = n_{i} mod n_{i+1}. In other words, n_{i+2} is the remainder when n_{i} is divided by n_{i+1}.

In Maple:

euc := proc(x::integer, y:: integer)
  if x < y then euc(y,x)                # swap
  elif y = 0 then x                     # terminate
  else euc(y, x mod y) fi               # reduce
end:

Theorem 1 Suppose x and y are integers, not both zero. Let M = {ax + by : a, b : Z}. Then GCD(x, y) = Minimum{g : g ∈ M, g > 0}.

Proof. The proof uses the fact that 0 ≤ x mod y <| y |. □

There is also an extended Euclidean algorithm which finds a and b so that ax + by = GCD(x, y). Try to write this in Maple. It should return the sequence a, b, GCD(x, y).

3.1 efficiency

Complexity of Euclid GCD over Z is O(n), where n is the size of the input, i.e. length in digits.

To show this, consider remainder sequence of length k, and pick x and y as small as possible with such a length. In order to do this, take all the quotients q_{i} to be 1. We get a backwards Fibonacci sequence as the worst case.

f_{i}, the i th number in the Fibonacci sequence has a remainder sequence of length i − 1 when divided by f_{i−1}, and no smaller number has a remainder sequence so long. Therefore the number of steps taken by the Euclidean algorithm is bounded by a multiple of the least i so that f_{i} is larger than max(| x |, | y |).

We have f_{i} < 2f_{i−1} for i bigger than 2, and so f_{i+1} > 3/2f_{i} > C(3/2)^{i}, for some C, and for i bigger than 2. So i < O(log(f_{i})), which gives us our result that the length of the remainder sequence starting with x and y is bounded by O(n), where n is the size of (x, y) measured in digits.

It is possible to get a better estimate for the Fibonacci sequence, using standard techniques for solving recurrence relations. It can be shown that f_{i} = c_{1}ϕ^{i} + c_{2}(−ϕ)^{−i}, where ϕ = (1 + 5^{1/2})/2, and c_{1}, c_{2} are 5^{−1/2}, −5^{−1/2} respectively.

4 Continued Fractions

See H. Davenport, Higher Arithmetic, (512.41 DAV).

Definition. 1 CF (q0, q1, ..., qn), the continued fraction with partial quotients q0, q1, ..., qn, is

  q0 + 1
      ----------------
      q1 + 1
           ------------
            ...
                  +1
                  ---
                   qn

CF (q_{0}, ..., q_{n}) = P_{n}(q_{0}, ..., q_{n})/Q_{n}(q_{0}, ..., q_{n}) where P_{n} and Q_{n} are polynomials in Z[q_{0}, ..., q_{n}]. CF (q_{0}, ..., q_{n}) is called the n th convergent of the sequence of partial quotients q_{0}, ..., q_{n}, ....

floor(x) is the largest integer n with n ≤ x.

Definition. 2 Suppose α : R. Define cf(α) = floor(α) if α = floor(α), and otherwise let cf(α) = floor(α) + 1/(cf(1/(α − floor(α))).

This defines a sequence of partial quotients q_{0}, .... The sequence terminates if α is rational and otherwise it goes on forever. The convergents are good approximations to α.

Lemma 1 Q_{n}(q_{0}, ..., q_{n}) = P_{n−1}(q_{1}, ..., q_{n})

P_{n}(q_{0}, ..., q_{n}) = q_{0}P_{n−1}(q_{1}, ..., q_{n}) + P_{n−2}(q_{2}, ..., q_{n})

Proof. CF (q_{0}, ..., q_{n}) = P_{n}/Q_{n} = q_{0} + 1/CF (q_{1}, ..., q_{n}) = q_{0} + Q_{n−1}(q_{1}, ..., q_{n})/P_{n−1}(q_{1}, ..., q_{n}). Put over common denominator.

□

Theorem 2 P_{n+1}(q_{0}, ..., q_{n+1}) = q_{n+1}P_{n}(q_{0}, ..., q_{n}) + P_{n−1}(q_{0}, ..., q_{n−1})

Q_{n+1}(q_{0}, ..., q_{n+1}) = q_{n+1}Q_{n}(q_{0}, ..., q_{n}) + Q_{n−1}(q_{0}, ..., q_{n−1})

Proof. By induction on n.

□

Theorem 3 For any real α there exists p/q in Q so that | α − p/q |< 1/q^{2}

Proof. Clear if α is rational. Suppose not. Define x mod 1 to be x − floor(x). Consider

M_{n} = {qα mod 1 : 0 < q ≤ n} ∪ {0, 1}

M_{n} consists of n + 2 distinct numbers, all in [0, 1]. By the pigeon hole principle, there are two distinct numbers x and y in M_{n} which differ by no more than 1/n.

| q_{1}α − p_{1} − (q_{2}α − p_{2}) |≤ 1/n

So | α − (p_{2} − p_{1})/(q_{1} − q_{2}) |≤ 1/(n | q_{1} − q_{2} |) < 1/(q_{1} − q_{2})^{2}

□

The proof shows something stronger: infinitely many of these good approximations. How to find them?

Theorem 4 | P_{n}/Q_{n} − α |= Minimum{| p/q − α :| q |≤ Q_{n} | | P_{n}/Q_{n} − α |< 1/(Q_{n}Q_{n+1}) < 1/Q_{n}^{2} where P_{n}/Q_{n} is n the convergent to α.

Proof omitted. See Davenport.

This theorem means that the n convergents are optimal rational approximations to α.

4.1 Coursework 1: M70

Write a procedure which computes the partial quotients for the continued fraction representation for given real α. Also give a procedure to compute the nth convergent to α. Consider α = (n)^{1/m}, for n and m in the range given to you by dsr. Find the first 51 partial quotients, q_{0}..q_{50} for each such number. Find the best available rational approximation for each α, and tabulate the differences between each α and its best approximation. Marking: Simplicity and apparent correctness of code 50; Explanation 25; Tests and results 25.

Due Friday March 9. Hand in to Mathematics Departmental Office before 3.30 pm.

5 Diophantine Equations

A Diophantine problem is one in which we are looking for whole number solutions to some condition.

5.1 Pell Equation

x^{2} − Dy^{2} = 1

5.2 Fermat’s last theorem and the ABC conjecture

x^{n} + y^{n} = z^{n}, n > 2

Conjecture 1 ABC conjecture. There exists a real number κ so that if A, B, C are natural numbers which are relatively prime and such that A + B = C, then

C ≤ (∏_{p|ABC,p:prime} p)^{κ}

6 LLL and PSLQ algorithms

6.1 Gram Schmidt

Gram Schmidt orthogonalisation process Let x^{T} mean the transpose of vector or matrix x. Recall that two non zero vectors a and b in R^{n} are orthogonal if and only if the vector product (i.e.

dot product, scalar product) a ∘ b = ab^{T} is zero. For example, (1, 1) ∘ (1, −1) = 0. Orthogonality is a nice property of a basis of a vector space.

Let v_{1}, ..., v_{k} be a list of vectors in R^{n} which is a basis for a vector space S(v_{1}, ..., v_{k}). The Gram Schmidt process produces an orthogonal basis v_{1}∗, ..., v_{k}∗ for this vector space as follows.

Define v_{1}∗ := v_{1}, and B_{1} := v_{1} ∗ ∘v_{1}∗, where ∘ means vector product. Now continue recursively, for all i > 1, for all j < i set μ_{ij} := v_{i} ∘ v_{j} ∗ /B_{j} (μ_{ij} is the projection of v_{i} in direction v_{j}∗, scaled by v_{j}∗.) Then v_{i}∗ := v_{i} − ∑_{j<i} μ_{ij}v_{j}∗, and B_{i} := v_{i} ∗ ∘v_{i}∗. So v_{i}∗ is what is left when we take away whatever can be expressed in terms of v_{1}∗, ..., v_{i−1}∗. Note that

the output (v_{1}∗, ..., v_{k}∗) need not be unit vectors.

6.2 Correction to proof of theorem 3

Add 0 and 1 to M_{n}.

M_{n} = {qα mod 1 : 0 < q ≤ n} ∪ {0, 1}

Then M_{n} consists of n + 2 numbers all in [0, 1].

6.3 More about the coursework 1

You will be finding continued fractions for n^{1/m}. Get your number n from the class list pinned to door of 1w2.2a. Pick bound M. If you are not on the list, pick a new number less than 50. You should hand in

1. Code for procedures which compute continued fraction for n quotients for real α, and kth convergent for same.

2. Some tests of the above, attempting to demonstrate correctness

3. For your n and for m from 2 to M, find partial quotients q_{0}, ..., q_{50}. Find largest q_{k}, and report k − 1 convergent, which should be the best approximation for n^{1/m}, and the number of decimal places of accuracy of the approximation.

4. What precision do you need? Discuss and explain.

5. Decide how to measure goodness of approximation. (See Thue Siegel Roth theorem.) Give best overall approximation which you found.

6.4 Quick review of linear algebra

Let v_{1}, ..., v_{k} be vectors in R^{n}.

Definition. 3 v_{1}, ..., v_{n} are linearly independent if

c_{1}v_{1} + ... + c_{n}v_{n} = 0

is impossible with c_{1}, ..., c_{n} real numbers which are not all zero.

Definition. 4 V ⊆ R^{n} is a vector space if V is closed under vector addition and subtraction, and V is closed under multiplication by real scalars.

Example: plane through origin.

Definition. 5 Span(v_{1}, ..., v_{n}) is the set of vectors which can be written in the form c_{1}v_{1} + ... + c_{n}v_{n} where c_{1}, ..., c_{n} are real numbers.

Definition. 6 (v_{1}, ..., v_{k}) is a basis for a vector space V if V = Span(v_{1}, ..., v_{k}) and v_{1}, ..., v_{k} are linearly independent.

Notation: If v_{1} = (x_{1}, ..., x_{n}) and v_{2} = (y_{1}, ..., y_{n}) then v_{1} ∘ v_{2} = x_{1}y_{1} + ... + x_{n}y_{n}. | v |= (v ∘ v)^{1/2}. This is Euclidean length.

6.5 Lattice Reduction

Try ?lattice; in maple.

Advice: search for Lenstra AND LLL in yahoo, and read around the subject. Let (b_{1}, ..., b_{k}) be a list of linearly independent vectors in R^{k}. A lattice in R^{k} is a subset L(b_{1}, ..., b_{k}) of the form {z_{1}b_{1} + ... + z_{k}b_{k}}, with z_{1}, ..., z_{k} are integers. We will

say that (b_{1}, ..., b_{k}) is a basis of L(b_{1}, ..., b_{k}).

The problem of finding a shortest non zero vector in a given lattice is considered to be very difficult. There is, however, a good method, the LLL algorithm, for finding an approximation to a shortest vector. This is based on the Gram Schmidt method, which was discussed earlier. The essential geometric idea is that if the basis were orthogonal, the smallest vector in it would also be the smallest in the lattice. So the algorithm follows the Gram Schmidt methods as closely as possible, but using integral coefficients. (We can not subtract arbitrary multiples of one vector from another since that might take us outside the lattice.)

The following definition (due to Lenstra, Lenstra, Lovasz) uses the notation of the Gram Schmidt method.

Definition. 7 A basis (b_{1}, ..., b_{k}) is LLL reduced if it satisfies the conditions:

1. | μ_{ij} |≤ 1/2 for all 1 ≤ j ≤ i ≤ k

2. | b_{i}^{∗} + μ_{ii−1}b_{i−1}^{∗} |^{2}≥ (3/4) | b_{i−1}^{∗} |^{2}, for 1 < i ≤ k.

It is intended to capture the idea that the vectors are approximately orthogonal.

Theorem 5 Suppose basis (v_{1}, ..., v_{n}) has Gram Schmidt basis (v_{1}^{∗}, ..., v_{n}^{∗}). Let m be any non zero vector in the lattice generated by (v_{1}, ..., v_{n}). Then | m |≥ Minimum(| v_{i}^{∗} |: 1 ≤ i ≤ n)

Proof. m = a_{1}v_{1} + ... + a_{n}v_{n}, where a_{1}, ..., a_{n} are integers, not all zero. Also m = c_{1}v_{1}^{∗} + ... + c_{n}v_{n}^{∗}. Pick j maximal so that a_{j} ≠ 0.

(v_{1}, ..., v_{j}) Gram Schmidt’s into (v_{1}^{∗}, ..., v_{j}^{∗}). Thus a_{1}v_{1} + ... + a_{j}v_{j} = m = c_{1}v_{1}^{∗} + ...c_{j}v_{j}^{∗}. We also have v_{j} = v_{j}^{∗} + ∑_{k<j} μ_{jk}v_{k}^{∗} (v_{1}, ..., v_{j−1}) Gram Schmidt’s into (v_{1}^{∗}, ..., v_{j−1}^{∗}). Thus m = a_{j}v_{j}^{∗} + T erms in (v_{1}^{∗}, ..., v_{j−1}^{∗}) = c_{j}v_{j}^{∗} + T erms in (v_{1}^{∗}, ..., v_{j−1}^{∗}), and since the basis vectors are

linearly independent, a_{j} = c_{j}. Thus c_{j} is a non zero integer, and | m |≥| v_{j}^{∗} |.

□ Another proof can be given via taking the dot product of m and v_{j}^{∗}. The result we need is:

Theorem 6 Let L be a lattice and let (b_{1}, ..., b_{k}) be a reduced basis of L. Let x be any non zero vector in L. Then | x |≥ 2^{−(n−1)/2} | b_{1} |.

Proof.

Perform Gram Schmidt on (b_{1}, ..., b_{k}) to get (b_{1}^{∗}, ..., b_{k}^{∗}). We know (by the previous theorem) that the lengths of the orthogonal vectors give us a lower bound on the length of any vector in the lattice. Since the basis is reduced, we also have

| b_{i}^{∗} + μ_{ii−1}b_{i−1}^{∗} |^{2}≥ (3/4) | b_{i−1}^{∗} |^{2}, for 1 < i ≤ k.

Thus | b_{i}^{∗} |^{2}≥ (1/2) | b_{i−1}^{∗} |^{2}

So for each i

| b_{i}^{∗} |≥| b_{1} | 2^{(n−1)/2}, and therefore this is also a lower bound for length of non zero vectors in the lattice. □ Here is a version of LLL. LLL algorithm Input (v_{1}, ..., v_{n}) Apply Gram Schmidt to get (v_{1}^{∗}, ..., v_{n}^{∗}) and (μ_{ij}) (so that we can test if (v_{1}, ..., v_{n}) is reduced) While (not reduced (v_{1}, ..., v_{n}))

1. Reset each v_{i} to v_{i} − ∑_{j<i} nint(μ_{ij})v_{j}, for i = 2..., n

2. Do Gram Schmidt process again getting (v_{1}∗, ..., v_{n}∗), and (μ_{ij})

3. If | v_{i}^{∗} + μ_{ii−1}v_{i−1}∗ |^{2}< 3/4 | v_{i−1}^{∗} |^{2}, for some i, then swap v_{i} and v_{i−1} for first such i

4. Apply Gram Schmidt again, to get (v_{1}^{∗}, ..., v_{n}^{∗}) and (μ_{ij})

Note similarity with Euclid GCD. Essentially we take little things away from big things as far as possible, and then we swap, and continue until finished. Proof of termination (omitted here) is much harder for LLL than for GCD , but LLL does always finish in polynomial time. Complexity is O(n^{6}(log | v_{ij} |)^{3}. Much slower than GCD over Z, but still usable on medium size problems. Dependence on dimension n is worse than on size of coefficients, as we would expect. Many optimizations and variations of the above technique are possible.

6.6 Integer relation finding

See papers in:

www.farcaster.com/papers/sm-thesis/node5.html
(LLL exposition)
 www.nersc.gov/~dhb/
(Analysis of PSLQ, an integer relation finding algorithm)
www.cecm.sfu.ca/personal/borwein/
www.cecm.sfu.ca/Projects/IntegerRelations/index.html
(Examples in Maple)
www.computer.org/cise/cs2000/c1toc.htm
(Integer relation finding)

Definition. 8 Suppose x = (x_{1}, ..., x_{n}) is a vector in R^{n}. An integer relation for x is a vector m = (m_{1}, ..., m_{n}) of integers, not all zero, so that xm^{T} = 0.

An n by n matrix Q is unitary if QQ^{T} = I. Let U(n, R) denote the set of unitary n dimensional matrices with real coefficients.

Let GL(n, Z) denote the set of n by n matrices with integral coefficients, and with determinant 1. GL(n, Z) is a group.

A matrix is lower trapezoidal if it has more rows than columns and if it is all 0 above the diagonal.

Definition. 9 Suppose x = (x_{1}, ..., x_{n}) is a unit vector in R^{n}, and that no coordinate entry of x is 0. For 1 ≤ j ≤ n, define

s_{j}^{2} = ∑_{j≤k≤n} x_{k}^{2} Define the n by n − 1 lower trapezoidal matrix H_{x} = (h_{ij} by h_{ij} = 0 if 1 ≤ i < j ≤ n − 1 h_{ij} = s_{i+1}/s_{i} if 1 ≤ i = j ≤ n − 1 h_{ij} = −x_{i}x_{j}/(s_{j}s_{j+1}) if 1 ≤ j < i ≤ n

Lemma 2 1. H_{x}^{T} H_{x} = I_{n−1}, i.e. the columns of H_{x} are orthogonal, and unit.

2. | H_{x} |= \sqrt{(}n − 1)

3. xH_{x} = 0

The columns of H_{x} are a basis for x^{⊥}, the space orthogonal to x.

Lemma 3 Let P_{x} = H_{x}H_{x}^{T}. Then

1. P_{x}^{T} = P_{x}

2. P_{x} = I_{n} − x^{T} x

3. P_{x}^{2} = P_{x}

4. | P_{x} |= \sqrt{(}n − 1)

5. P_{x}z^{T} = z^{T} for any z in x^{⊥}

6. P_{x}m^{T} = m^{T} for any integer relation for x.

Theorem 7 Suppose that m is an integer relation for x. Let H = AH_{x}Q, where A is in GL(n, Z) and Q is in U(n − 1), and suppose that H is lower trapezoidal, and each diagonal element h_{j,j} ≠ 0. Then

| m |≥ Minimum(1/ | h_{ii})

Proof. P_{x} = H_{x}H_{x}^{T} = I_{n} − x^{T} x. P_{x}m^{T} = m^{T}, as in lemma above. AP_{x}m^{T} = Am^{T} AH_{x}H_{x}^{T} m^{T} = Am^{T}. We have AH_{x} = HQ^{T}, since QQ^{T} = I_{n−1}. So HQ^{T} H_{x}^{T} m^{T} = Am^{T}, and Am^{T} is a column of integers, not all 0. Let Q_{H,j} be the j th row of Q^{T} H_{x}^{T}. This is a unit vector since Q^{T} H_{x}^{T} H_{x}Q = I_{n−1}. Let A_{j} be the j th row of A. Pick j minimal with A_{j}m^{T} ≠ 0. Q_{H,k}m^{T} = 0 for k < j, since j was chosen minimal. A_{j}m^{T} = h_{jj}Q_{H,j}m^{T}. So | m |≥ 1/h_{jj}. □

The idea of this theorem is that H is a transformation of H_{x}, using A and Q. We get a lower bound on the size of any integer relation for x by looking at the diagonal elements of H. Any matrices A and Q can be used here, provided only that H is lower trapezoidal. As long as we choose the A and Q so that the lower bound is increasing, we are making progress.

PSLQ algorithm

Pick γ > 2/(3)^{1/2}.

Set H = H_{x}

While (no element of x or diagonal of H is zero)

1. Hermite reduce H, by integral row operations, (subtracting integral multiples of the i th row from rows below i), so that in each column the diagonal has at least twice absolute value of the other elements. H := AH, x := xA^{−1}, where A in GL(n, Z) represents the reduction. The reduced matrix is still lower trapezoidal.

2. Find r so that γ^{r} | h_{rr} | is maximal. Interchange row r and row r + 1 of H. Let R represent this interchange. H := RH, x := xR^{−1}. Let α = H(r, r), β = H(r+1, r), λ = H(r+1, r+1), δ = (β^{2}+λ^{2})^{1/2}

3. H may not be lower trapezoidal at this point, since the interchange may have created a corner above the diagonal. Do column operations, represented by unitary matrix Q to get rid of this corner. H := HQ. Q is in U(n − 1), and may be defined by the following: Q(r, r) = β/δ, Q(r, r + 1) = −λ/δ, Q(r + 1, r) = λ/δ, Q(r + 1, r + 1) = β/δ, and otherwise Q is the n − 1 by n − 1 identity.

At the beginning and end of each loop, H is lower trapezoidal. Termination implies that an integer relation has been found. Non termination, means that any relation must satisfy | m |≥ Min(1/(| h_{ii} |), and this number increases at each loop. Note similarity with Euclidean algorithm.

An implementation in maple , written by Simon Langley, can be found (together with some examples):

~masdr/public_html/pslq.ms
~masdr/public_html/pslq.mws

6.7 A sample application of PSLQ

Suppose we have the ability to approximate a real number α with high accuracy, and we believe α is a root of a polynomial with integral coefficients of degree n.

a_{n}α^{n} + … a_{0} = 0

with (a_{n}, …, a_{0}) integers. Then (a_{n}, …, a_{0}) is an integer relation for x = (α^{n}, …, α, 1). So we can use PSLQ to search for the coefficients.

On the other hand if we believe that α is not algebraic of degree n, we can use PSLQ to exclude the possibility of polynomials with coefficients which are all below some bound, M. To do this, just do PSLQ until the bound produced by the theorem above exceeds nM.

6.8 Programming exercises

These will be discussed in the problems classes. The object is to learn how to use matrices in maple, and, in general, the linear algebra package. An important point: if B is a matrix, A:=B sets address of A to address of B, so that afterwards A and B are not independent. If you want an independent copy, you need to declare A and then write a loop or so to copy B into A.

1. Write Gram Schmidt procedure. (There is one already in the linear algebra package. To see the code for this and other procedures in maple, use interface(verboseproc=2), print(procedurename). Refer to help.)

2. Write procedure which computes H_{x} from x for pslq.

7 Computation in Z(x_{1}, …, x_{n})

7.1 Some Algebra

7.2 Introduction

Definition. 10 Suppose K is a set. A polynomial in x with coefficients in K is an expression of the form a_{d}x^{d} + a_{d−1}x^{d−1} + ... + a_{0} where d is a natural number, and a_{d}, ..., a_{0} are all elements of K.

Notation: K[x] is the set of polynomials in x with coefficients in K.

Definition. 11 A power series in x with coefficients in K is an expression of the form ∑_{n=0}^{∞} a_{n}x^{n} = a_{0} + a_{1}x + a_{2}x^{2} + ... where a_{0}, a_{1}, ... are all in K.

Notation: K[[x]] is the set of power series in x with coefficients in K. Normally it is assumed that K is a type with +, −, ∗ defined, and in that case we can define +, −, ∗

uniformly in K[x] from the definitions in K. We do, however, need to make some assumptions about how +, −, ∗ behave.

7.3 Basic Definitions

Definition. 12 In the following a ring is a structure (D, +, ∗, 0, 1, =) where

• D is a collection of objects

• +, ∗ are operations from D × D to D

• 0 and 1 are in D, and 0 ≠ 1

• = is the usual equality on D

• (D, +, 0) is an Abelian group

• ∗ is associative and commutative, and 1 ∗ x = x for all x in D

• x ∗ (y + z) = x ∗ y + x ∗ z for all x, y, z in D

The full algebraic name of this idea is “commutative ring with unit”, but we will just call it a ring in the following.

Definition. 13 An integral domain is a ring in which x ∗ y = 0 only happens when either x = 0 or y = 0.

Definition. 14 A field is a commutative ring with unit in which every non zero object has a multiplicative inverse.

Every integral domain extends uniquely to a field of fractions. For example, Q is the field of fractions of Z.

If I is an integral domain, I[X] is the integral domain of polynomials in X with coefficients in I. I(X) is the field of fractions. Similarly, K[x_{1}, ..., x_{n}] is multivariate polynomials in x_{1}, ..., x_{n} with coefficients in K; and K(x_{1}, ..., x_{n}) is the field of fractions.

Definition. 15 x is a unit in ring D if there is a y in D so that x ∗ y = 1.

Definition. 16 Elements x and y in a ring are associates if there is a unit u in the ring so that x ∗ u = y

Definition. 17 We will say that x is a factor of y in a ring D, and write this as x | y if there is a z in D so that x ∗ z = y.

Notation: ∧ means logical AND. ∨ means logical OR. ¬ means logical negation. → means logical implication. ↔ means if and only if. ∀ means for all. ∃ means there exists.

Definition. 18 Suppose K is an integral domain and a and b are in K. An object g in K is a gcd of a and b if

g | a ∧ g | b, and

(∀h ∈ D)((h | a ∧ h | b) → h | g)

In other words if g is a gcd of a and b, then g divides both a and b, and g is as large as it can be with this property in the sense that g is divisible by anything in K which divides both a and b.

According to this definition both -4 and 4 would be gcds of 12 and 16 over Z. When K has an ordering, we usually define the GCD to be the largest object which satisfies the gcd definition given above.

Theorem 8 If g_{1} and g_{2} are gcds of x and y in integral domain K, then g_{1} and g_{2} are associates.

A little later in the course we will look at how to compute polynomial gcds. At this point, we do not even have a proof that they exist.

Definition. 19 Let K be an integral domain in which gcds exist. For p in K[x], define content(p) = gcd({Coefficients(p)}). Define P rimitiveP art(p) = p/content(p).

A polynomial is primitive if its content is 1.

Definition. 20 A non unit polynomial p is irreducible in K[x] if p = q ∗ r is impossible for q and r in K[x], unless either q or r is a unit.

Theorem 9 Z[x_{1}, ..., x_{n}] has unique factorisation.

(Every non zero non unit polynomial in this domain can be expressed as a product of irreducible factors in this domain. Furthermore the product is unique except for the order of the factors, and multiplication by units.)

This theorem implies that gcds are defined in Z[x_{1}, ..., x_{n}].

7.4 Ring morphisms

Definition. 21 Suppose K_{1} = (D, +, ∗, zero(K_{1}), one(K_{1}), =) and K_{2} = (E, ⊕, ⊗, zero(K_{2}), one(K_{2}), =) are rings. A ring morphism from K_{1} to K_{2} is a map

f : D → E

so that

1. f(x + y) = f(x) ⊕ f(y)

2. f(x ∗ y) = f(x) ⊗ f(y)

3. f(one(K_{1})) = one(K_{2}).

4. f(zero(K_{1})) = zero(K_{2})

Comments: In the first item above, the first addition x + y is going on in K_{1}, and the second addition f(x) ⊕ f(y) is going on in K_{2}. So this says that if you first add in K_{1} and then apply f, you will get the same result as if you first applied f to x and y and then added the results in K_{2}.

The fourth condition is a consequence of the first three Ring morphisms are also called ring homomorphisms. Ring morphisms are to rings what group homo-

morphisms are to groups.

The category of rings consists of all rings and all ring morphisms.

Example 3 The natural map f : Z → Z_{p}, where f(n) is n modulo p, is a ring morphism.

7.5 Representation of Rings

A representation of a ring is a collection of expressions which stand for objects in the ring. A representation is normal if 0 has only one representation. A representation is canonical if every object in the ring has only one representation.

If E_{1} and E_{2} are two expressions in a representation, we will say that E_{1} = E_{2} if E_{1} and E_{2} stand for the same object. For example, 0 + 0 = 0, even though 0 + 0 and 0 are not the same as expressions.

We will say that a representation has decidable identity problem if there is an algorithm to decide whether or not E_{1} = E_{2}, given E_{1} and E_{2}.

7.6 Integers of arbitrary size

Choose a base b, and represent integers, canonically, as either 0 or as a string d_{k}d_{k−1}...d_{0} or as a signed string −d_{k}d_{k−1}...d_{0} with 0 ≤ d_{i} < b.

Theorem 10 Division Theorem for Z.

For any a and b ≠ 0 in Z, there exists unique q and r in Z so that a = q ∗ b + r∧ | r |<| b | ∧ | q ∗ b |≤| a |

We will write the number r asserted to exist by this theorem as rem(a, b).

Definition. 22 a and b are relatively prime if gcd(a, b) = 1.

7.7 Rational numbers

We can represent rational numbers as p/q where gcd(p, q) = 1 and q > 0. This is a canonical representation.

Addition, subtraction, multiplication and division of rational numbers requires computation of gcds in Z, in order to put the result into canonical form.

7.8 Representation of Polynomials

7.8.1 Polynomials in K[x]

Suppose we have defined a canonical form for ring K. We can write a polynomial p(x), of degree d, in K[x], as an expression

p(x) = a_{0}x^{d} + ... + a_{1}x + a_{d} with the coefficients a_{0}, ..., a_{d} in canonical form, and a_{0} ≠ 0. On the data structure level, we may then either represent p(x) as an array of length d+1 of the coefficients (a_{0}, a_{1}, ..., a_{d}) or as a list of pairs of coefficients and powers, for the non zero coefficients, i.e.

{(a_{0}, d), ...}. The array representation is called dense and the list representation is called sparse.

7.9 Multivariate Polynomials

Suppose we are in K[x_{1}, ..., x_{n}]. There are a lot of choices as to how to write these polynomials.

7.9.1 Recursive sparse

The recursive sparse representation does this as K[x_{1}, ..., x_{n−1}][x_{n}]. That is, we consider a polynomial to be in x_{n} with coefficients in K[x_{1}, ..., x_{n−1}]. So p in K[x_{1}, ..., x_{n}] is written as

p = a_{0}x_{n}^{d} + ... + a_{d}

where the coefficients a_{0}, ..., a_{d} are all in canonical form in K[x_{1}, ..., x_{n−1}], and a_{0} is not 0, unless p = 0. This gives us a canonical form, which depends upon initially ordering the variables by importance x_{1} ≺ x_{2} ≺ ... ≺ x_{n} Of course, if we change the order of importance of the variables, the canonical forms change. This representation is also called lexicographic recursive. If y is the most important variable which actually occurs in polynomial p, we can write p as a_{0}y^{d} + a_{1}y^{d−1}... + a_{d} with a_{0} ≠ 0. In this case we will define lv(p) = y lc(p) = a_{0} deg(p) = d reduct(p) = a_{1}y^{d−1} + ... + a_{d} Assume deg(p) > 0. p could then be represented as the list {{{lv(p), deg(p)}, lc(p)}, reduct(p)}. So, in terms of Lisp, the car of a polynomial is the leading term and and cdr is the reduct.

7.9.2 Other canonical forms in Q[x_{1}, ..., x_{n}]

Definition. 23 A power product is a term obtained by multiplying together the variables to various non negative powers.

x_{1}^{d_{1}}...x_{n}^{d_{n}} with d_{i} ≥ 0 for i = 1, ..., n. The sum of the powers, d_{1} + d_{2} + ... + d_{n}, is the total degree of a power

product.

Definition. 24 A monomial of K[x_{1}, ..., x_{n}] is the product of a power product and an element of the ring K

Definition. 25 Suppose we order the variables x_{1} ≺ x_{2} ≺ ... ≺ x_{n}. Let a = c_{1}x_{1}^{a_{1}}x_{2}^{a_{2}}...x_{n}^{a_{n}} and b = c_{2}x_{1}^{b_{1}}x_{2}^{b_{2}}...x_{n}^{b_{n}} be monomials. Then we will say a ≺ b if a_{n} < b_{n} or if a_{n} = b_{n} but a_{n−1} < b_{n−1} or if a_{n} = b_{n} and... a_{1} = b_{1} but a_{1} < b_{1}. This is called lexicographic ordering on the monomials.

Each ordering of the variables induces a different lexicographic ordering.

Definition. 26 Another common ordering on monomials is by total degree, and then, in case of a tie, lexicographically.

We can get various canonical forms by writing our polynomials as sums of monomials, in decreasing order, depending on what order we adopt. This can be represented as lists of representations of monomials. This is called a distributed representation.

A common form is obtained by using the total degree and then lexicographic order on the monomials. Another common form is to use total degree and then reverse lex.

7.9.3 Fractions of polynomials

We can represent fractions of polynomials as ratios of polynomials with gcd 1.

In an integral domain, any two gcds of p and q are associates. We will normally talk about the gcd in this case, not distinguishing between associates.

7.10 Power Series

In terms of expressions, we can write a power series as ∑_{n=0}^{∞} a_{n}x^{n} where a_{n} is some expression depending on n. A truncated power series is an expression of the form ∑_{n=0}^{k} a_{n} where k is some natural number. One truncated power series ∑_{n=0}^{k} a_{n}x^{n} extends another power series ∑_{n=0}^{j} b_{n}x^{n} if j < k and a_{n} = b_{n} for all n ≤ j. Of course we do not literally have to write down all infinitely many terms of a power series, as long as we

have a rule for generating the coefficients. It would seem that we could construct a truncation of the power series, which we could use as an approximation, and keep some process ready to extend the truncation if needed.

Definition. 27 We can represent an infinite power series as (T, P) where T is a truncated power series, and P is a process which produces another infinite power series representation (T ∗, P ∗) where T ∗ extends T.

The idea here is that we can see any finite part, as much as we like, of the power series. If we want more of it there is a process in the background which will get it. These are sometimes called lazy series.

7.11 Polynomial GCD

7.12 Gauss lemma and GCDs in Z[x]

The relationship between factorisation (and GCDs) in Z[x] and Q[x] was studied by Gauss. Gauss defined the content of a polynomial to be the GCD of its coefficients. A polynomial is said to be primitive if its content is 1. If A is a polynomial, pp(A) = A/content(A). pp(A) is called the primitive part of A. (Note that these definitions generalise to K[x] where K is any ring in which there exist GCDs.)

Gauss Lemma For A, B in Z[x], if A and B is primitive, so is AB. Proof. By contradiction. Suppose A and B are primitive and AB is not primitive. So content(AB) = c,

where c is not a unit. Let k be a prime so that k | c. Let f_{k} : Z[x] → Z_{k}[x] be the ring morphism which reduces coefficients mod k. That is

f_{k}(a_{d}x^{d} + ... + A_{0}) = (a_{d}modk)x^{d} + ...a_{0}modk Since k does not divide content(A) and does not divide content(B), we have f_{k}(A) ≠ 0, f_{k}(B) ≠ 0. But Z_{k}[x] is an integral domain, f_{k}(AB) = 0 since k does divide content(AB); and also f_{k}(AB) =

f_{k}(A) ∗ f_{k}(B) since f_{k} is a ring morphism. So we have a contradiction.

□

Suppose A and B are in Z[x].

Corollary 1 pp(AB) = pp(A)pp(B)

Corollary 2 content(AB) = content(A)content(B)

Corollary 3 A | B iff (content(A) | content(B) and pp(A) | pp(B))

Corollary 4 GCD(A, B) = GCD(content(A), content(B)) ∗ GCD(pp(A), pp(B))

Corollary 5 if A has a non trivial factorisation in Q[x], then A also has a non trivial factorisation in Z[x].

7.13 Gauss lemma and GCDs in Z[x_{1}, ..., x_{n}]

We say that an integral domain has unique factorisation if, as in Z, every non zero non unit object can be expressed as a product of irreducible factors; and if such factorisation is unique, except for change of order of factors and multiplication by units.

We will use the following theorem without proof.

Theorem 11 Z[x_{1}, ..., x_{n}] has unique factorisation.

The Gauss lemma can be generalised (with only small change to the proof) from Z[x] to K[x], where K is an integral domain with unique factorisation. In particular, consider K[x_{n}], where K is Z[x_{1}, ..., x_{n−1}].

Gauss Lemma 2 For A, B in Z[x_{1}, ..., x_{n−1}][x_{n}], if A and B are primitive, so is AB.

All the corollaries also generalise.

7.14 Problems

1. Prove that if k is a prime integer then Z_{k} is an integral domain.

2. Prove that if I is an integral domain, so is I[x].

3. Write (r^{2} + s)(r − s + 1) in Z[r][s].

4. Represent (r^{2} + s)(r − s + 1) in sparse recursive form in Z[r][s].

5. Prove that if a, b, c are in Z[x] and ab = c, then con(a) ∗ con(b) = con(c), and pp(a) ∗ pp(b) = pp(c). Use the Gauss lemma, and show where it is used.

6. Show that for a, b in Z[x], GCD(a, b) = GCD_{Z}(con(a), con(b)) ∗ GCD_{Z[x]}(pp(a), pp(b)) Use previous problem.

7.15 GCDs in Z[x_{1}, ..., x_{n}]

We will look at a number of methods for finding these GCDs. This is a basic problem since it appears to be needed for addition and multiplication in Z(x_{1}, ..., x_{n}).

All of these methods follow the same recursive pattern, using the GCD corollary of the previous section. Suppose we have a way to find GCDs in Z[x_{1}, ..., x_{n−1}]. Let K = Z[x_{1}, ..., x_{n−1}]. Suppose given A and B in K[x_{n}]. Find content(A), content(B), pp(A), pp(B), using our existing method for GCDs in K. Then

GCD(A, B) = GCD(content(A), content(B)) ∗ GCD(pp(A), pp(B)).

We can find GCD(content(A), content(B)) by our existing method. We are left with finding GCD for pp(A) and pp(B). So we need a way to get GCDs for primitive polynomials. There are various ways of doing this.

wolog, in all the following, we assume A and B are primitive.

Assume in the following that K[x] is an integral domain.

7.15.1 Method 1

Given A and B primitive in K[x_{n}], extend K to field of fractions F r(K). Find the GCD of A and B in F R(K)[x_{n}] using the Euclidean algorithm. Suppose this GCD is G. The coefficients of G may be fractions of polynomials in K. Multiply by the least common multiple of the denominators, to clear all denominators, obtaining a polynomials g in K[x]. Then take pp(g). This is GCD(A, B) in K[x_{n}].

Example 4 A(x) = x^{8} + x^{6} − 3x^{4} − 3x^{3} + 8x^{2} + 2x − 5 B(x) = 3x^{6} + 5x^{4} − 4x^{2} − 9x + 21

A problem is that there may be many recursive GCD computations, and the coefficients may become huge. So we are led to consider the possibility of working directly in K, and not using fractions at all.

7.15.2 Pseudodivision technique

Finding remainders depends on division. So we define something which is like a remainder, but which will not need division which goes outside K.

Theorem 12 Suppose p(x) and q(x) are in K[x], and p(x) has degree m and q(x) has degree n. There exist c(x) and r(x) in K[x] with the degree of r(x) less than the degree of q(x), or r(x) = 0, so that

lc(q)^{m−n+1} ∗ p(x) = c(x) ∗ q(x) + r(x)

Definition. 28 The polynomial r(x) asserted to exist above will be denoted prem(p, q), prem(p, q) will be called the pseudo remainder of p with respect to q.

Method 2 for calculating gcds in K[x] is to use pseudoremainders instead of remainders in Euclid’s algorithm.

———————-

Definition. 29 The Pseudo Euclidean Algorithm for GCDs Assume A and B are primitive in K[x]. GCD(A, B) is defined by:

while(B ≠ 0)(R := prem(A, B), A := B, B := R) Then GCD = pp(A).

———————–

A problem with this is that the contents get huge. Try it on the previous example in Z[x]. Especially annoying is the fact that we do not want the contents. We are trying to get the GCD of primitive polynomials, and so we will end up with a primitive GCD.

Method 3 is to remove the contents at each stage. This has lots of recursive GCD computations.

7.15.3 Polynomial remainder sequences

We need some notation. We assume that A(x) and B(x) are in K[x] and that the degree of A(x) is at least as large as the degree of B(x). Define a polynomial remainder sequence for A and B to be a sequence of polynomials R_{0}, R_{1}, R_{2}, ..., R_{j} in K[x] satisfying:

1. R_{0} = A, R_{1} = B

2. For some α_{i}, Q_{i}, β_{i} in K[x],

α_{i}R_{i−1} = Q_{i}R_{i} + β_{i}R_{i+1}

3. prem(R_{j−1}, R_{j}) = 0

For the Euclidean polynomial remainder sequence using pseudo-remainders, α_{i} = r_{i}^{δ_{i}+1}, r_{i} = lc(R_{i}), δ_{i} = deg(R_{i−1}) − deg(R_{i}), β_{i} = 1.

For the primitive polynomial remainder sequence, α_{i} = r_{i}^{δ_{i}+1}, β_{i} = content(prem(R_{i−1}, R_{i})). A good deal of work has been done to try to predict part of the content of the polynomials in the sequence,

without actually having to calculate all the recursive GCDs. If you try some experiments, you will see that the multipliers (the r_{i}) often turn up later on in the contents. In fact, it can be shown that in the Euclidean sequence, we can always find a factor of r_{i−1} in R_{i+1}. Method 4 is to divide this out, giving us a reduced polynomial remainder sequence:

α_{i} = r_{i}^{δ_{i}+1}, β_{1} = 1, β_{i} = α_{i−1} for 2 ≤ i ≤ j. Quite often more than one factor of this kind appears. Method 5, the trial divisor method, reduces the

content by as large a power of r_{i−1} as divides evenly.

The ultimate thing of this kind is the subresultant polynomial remainder sequence which looks like this:

α_{i} = r_{i}^{δ_{i}+1}, β_{1} = (−1)^{δ_{i}+1}, β_{i} = −r_{i−1}Ψ_{i}^{δ_{i}} for 2 ≤ i ≤ j, with Ψ_{i} defined by: Ψ_{1} = −1 Ψ_{i} = (−r_{i−1})^{δ_{i−1}}Ψ_{i−1}^{1−δ_{i−1}} for 2 ≤ i ≤ j.

8 Numerical Solving by Newton’s Method

9 Newton’s Method

In this section F = (f_{1}, ..., f_{n}) : C^{n} → C^{n}, all the partial derivatives of all orders are assumed to exist, x = (x_{1}, ..., x_{n}), and J(x) is defined to be the n by n matrix of partial derivatives of F (x), i.e.

J(x) = (∂f_{i}/∂x_{j}). J(x) is called the Jacobian of F (x). The Jacobian, J(x), is an n dimensional derivative. It specifies the “slope” of the best linear approxima-

tion to F (x) at a point.

F (x) = F (x_{0}) + J(x_{0})(x − x_{0}) + O(| x − x_{0} |^{2}). The linear function F (x_{0}) + J(x_{0})(x − x_{0}) is called the linearisation of the function F (x) at x_{0}. A solution point of n equations in n unknowns, F (x) = 0, is said to be non singular if the Jacobian is

non singular at that point.

Suppose we want to solve F (x) = 0 and we have an initial approximation x = a to a solution. Newton’s method constructs a sequence of approximations to a solution by linearising the function at each point, and solving the linearised equation. The recurrence obtained is: x_{0} = a, x_{i+1} satisfies

F (x_{i}) + J(x_{i})(x_{i+1} − x_{i}) = 0

If J(x_{i}) is non singular, then:

x_{i+1} = x_{i} − J^{−1}(x_{i})(F (x_{i}))

This works nicely if our original guess x_{0} is sufficiently near a non singular solution.

Here is a famous theorem.

Theorem 13 Inverse Function Theorem. If F (x) = y and the Jacobian is non singular at x, then F has an inverse in some neighbourhood of y.

And here is a famous conjecture.

Conjecture 2 The Jacobian Conjecture *****. In either R^{n} or C^{n}, if F is a polynomial map and the Jacobian of F is never singular, then F is 1-1 and onto.

and another famous theorem:

Theorem 14 If F (a) = 0 and the Jacobian is continuous and non singular at x = a, there is a neighbourhood of a in which Newton’s method converges to the solution a no matter what initial point is chosen in the neighbourhood.

9.1 Interval Arithmetic

We can define an arithmetic of intervals as follows.

Definition. 30 [a, b] + [c, d] = [a + c, b + d] [a, b] ∗ [c, d] = [m, M], where m and M are the min and max respectively of (ac, ad, bc, bd).

We can represent numbers, r, as [r, r].

So [a, b] − [c, d] should be [a, b] + [−1, −1] ∗ [c, d].

So far we have:

Lemma 4 {xopy : x ∈ [a, b], y ∈ [c, d]} ⊆ [a, b]op[c, d] with op being any of +, −, ∗.

The same inclusion property holds for polynomials. We would like to make these intervals smaller if possible, while keeping the inclusion. So, for example, it is better to define [a, b]^{3} as [a^{3}, b^{3}], rather than as [a, b] ∗ [a, b] ∗ [a, b] as above. It seems that the right way to define [a, b]^{2} depends on the signs of a and b...

Definition. 31 A Q-box in R^{n} is a closed coordinate aligned box whose vertices have coordinates which are either −∞ or ∞ or rational numbers, expressed as bounded precision floats.

9.2 Validated numerical methods

Interval arithmetic can be used in solving systems of equations. If f is a function of n variables with values in R^{m}, and box is a box in R^{n} we will write f(box) to be the box in R^{m} obtained by the interval technique discussed above. So we are guaranteed that the image of box under f is contained in f(box). This already gives us a powerful method to exclude the possibility that f(x) = 0 has a solution in box.

0 not in f(box) → (¬(∃x ∈ box)f(x) = 0

This simple idea can be used with recursive subdivision of the box to make an elegant solver. Suppose we have a condition small(box) so that if this condition holds and 0 ∈ f(box), then we would accept box as an approximate solution.

The following algorithm either returns one such box or it returns the empty set, in which case there is no solution in the box.

Algorithm solve(f = 0, box)

1. if 0 is not in f(box), then return Ø

2. if 0 is in f(box) and small(box), then return box

3. Otherwise, subdivide box into box1 and box2, splitting at the middle of a side of maximal length. Return solve(f = 0, box1) if this is non empty. Otherwise, return solve(f = 0, box2).

Suppose we measure the length of a vector, v, of complex numbers by taking the square root of the dot product of v and conjugate(v).

If A is a matrix of complex numbers, let ∥A∥ be the maximum of the lengths of the rows of A. We have

| Aw |≤ ∥A∥ | w | n

, where A is square n by n and w is a vector of length n, with complex entries. We can also use interval arithmetic to verify convergence of Newton’s method.

Theorem 15 Suppose box is a box in C^{n}. Let newbox be (x − J^{−1}(x)F (x))(box).

Suppose:

1. newbox ⊂ box

2. For x and y in box, ∥I − J(x)J^{−1}(y)∥ < 1/2n

Then the Newton sequence started at any x_{0} ∈ box converges to a solution in box of F (x) = 0.

Proof. The iteration x_{n+1} = x_{n} − J^{−1}(x_{n})F (x_{n}) stays inside of box, since newbox ⊂ box. Consider one step. Draw a straight line, L, from x_{n} to x_{n+1}, and consider what happens to F (x) on this line. F (x_{n+1}) = F (x_{n}) + ∫_{L} J(x)dx Parameterise x(t) = x_{n}(1 − t) + x_{n+1}t, for 0 ≤ t ≤ 1. dx = −J^{−1}(x_{n})F (x_{n})dt F (x_{n+1}) = F (x_{n}) + ∫_{0}^{1}(J(x_{n}) + (J(x(t)) − J(x_{n}))(−J^{−1}(x_{n})F (x_{n}))dt ≤ Max | I − J(x)J^{−1}(x_{n}) || F (x_{n}) |, for x on the line L. Thus, by the second condition of the premise of the theorem, | F (x_{n+1}) |< 1/2 | F (x_{n}) |. This implies that the norm of F (x_{n}) tends to zero as the iteration proceeds, and thus the difference

between x_{n+1} and x_{n} tends to zero. Thus the Newton sequence converges.

□ We can prove convergence of the Newton sequence by verifying the two conditions using interval arith-

metic.

There are a number of other tests to verify Newton convergence. In conclusion, we have a powerful technique for solving n equations in n unknowns, provided that we begin with a good approximation to a non singular solution.

9.3 Coursework 2: Due Friday May 11

Hand in to Mathematics Office, before 4 pm. Write and test a procedure for Newton’s method in complex n space. The user should give

1. A vector α_{0} in C^{n} or R^{n}

2. A vector F (x) = (f_{1}(x), …, f_{n}(x)), with f_{i}(x) = f_{i}(x_{1}, …, x_{n}) : C^{n} → C, and F (x) : C^{n} → C^{n}.

3. A precision k. (Any number ϵ, with | ϵ |< 10^{−k} will be regarded as insignificantly small.)

4. A maximum number of iterations

The procedure should return either an α ∈ C^{n} so that | F (α) |< 10^{−k}, or a helpful error message of some kind.

Marking Scheme:

1. Design of procedure (20)

2. Simplicity and correctness of code (50)

3. Testing (30)

Here are some examples which can be used for testing, with various starting points.

• x^{3} − 1 : C → C. Try randomly chosen initial points in unit circle around 0 in C.

• (x_{1}^{3} − 1, x_{2}^{3} − x_{1}, ..., x_{n}^{3} − x_{n−1}) : C^{n} → C^{n}

• (e^{xz} − 2 ∗ y ∗ w, e^{xw} − 3, e^{yz} − 7, e^{yw} − 11) : C^{4} → C^{4}

• x^{4} + x^{3} − x − 1 : C → C.

10 Factorisation inZ and in Z[x]

10.1 The Tau conjecture

A computation of a number m in Z is a sequence x_{0}, …, x_{l} with x_{0} = 1 and x_{l} = m, and, for each i, x_{i} is obtained from two previous elements in the sequence

by either addition, subtraction or multiplication. Define τ(m) to be the length of the shortest computation which ends in m.

Conjecture 3 The τ conjecture. There is no constant C so that τ(k!) ≤ (log k)^{C} for k ∈ N

If the τ conjecture is false, i.e. if there is a quick way to compute k!, then there is also a quick way to factorise integers.

If n = pq with p < k < q, compute k! by quick computation x_{0}, ..., x_{l} and then get a computation of k! modn by x_{0} modn, ..., x_{l} modn = r and then take GCD(r, n) to get a factor of n. See Complexity and Real computation, Blum, Shub, Smale, Cucker, Springer, (510.095.9).

10.2 Polynomial factorisation

Height of polynomial. Landau Mignotte bound on height of factors. To factor p(x) in Z[x], find a root α, approximate with Newton’s method. Then use PSLQ to find a minimum polynomial q(x) defining α. Then q(x) is irreducible and q(x) divides p(x).